double distancia(double **pos, double lbox, int i, int j, double *dx, double *dy, double *dz)
{
	double ddx,ddy,ddz;
	double r2=0.;

	ddx=pos[i][0]-pos[j][0]; ddx-=dround(ddx/lbox)*lbox;
	ddy=pos[i][1]-pos[j][1]; ddx-=dround(ddy/lbox)*lbox;
	ddz=pos[i][2]-pos[j][2]; ddx-=dround(ddz/lbox)*lbox;
	r2=SQR(ddx)+SQR(ddy)+SQR(ddz);
	*dx=ddx;
	*dy=ddy;
	*dz=ddz;
	return r2;
}


double dipolarEnergyForces(double **pos, double lbox, double **force, int npart)
{
	double distancia(double **pos, double lbox, int i, int j, double *dx, double *dy, double *dz);
	double r,r2, r3, r4,dx,dy,dz, Co, Co2;
	double energia=0.0;
	double angPart1,angPart2;
	double fx,fy,fz;

	for(int i=0;i<npart;i++){
		for(int k=0;k<3;k++){
			force[i][k]=0.0;
		}
	}

	for(int i=0;i<npart-1;i++){
		for(int j=i+1;j<npart;j++){
			r2=distancia(pos, lbox, i, j, &dx, &dy, &dz);
			r=sqrt(r2);
			r3=r2*r;
			r4=r3*r;
			Co=dz/r;
			Co2=Co*Co;
			energia+=(1.-3.*Co2)/(r3);

			angPart1= (5.0*Co2-1.0)/r4;
			angPart2= (5.0*Co2-3.0)/r4;
			fx=angPart1*dx/r;
			fy=angPart1*dy/r;
			fz=angPart2*dz/r;

			force[i][0]+=fx;
			force[i][1]+=fy;
			force[i][2]+=fz;

			force[j][0]-=fx;
			force[j][1]-=fy;
			force[j][2]-=fz;
		}
	}
	return energia;
}

double dipolarEnergyForcesReplication(double **pos, double lbox, double **force, int npart)
{
	double distancia(double **pos, double lbox, int i, int j, double *dx, double *dy, double *dz);
	double r4,r3,r2, r,dx,dy,dz, Co, Co2;
	double dxn,dyn,dzn;
	double angPart1,angPart2;
	double fx,fy,fz;
	double energia=0.0;
	int ncut=30;
	int ncut2=ncut*ncut;
	int nq;

	for(int i=0;i<npart;i++){
		for(int k=0;k<3;k++){
			force[i][k]=0.0;
		}
	}

	for(int i=0; i<npart; i++){
		for(int j=0; j<npart; j++){

			dx=pos[i][0]-pos[j][0];
			dy=pos[i][1]-pos[j][1];
			dz=pos[i][2]-pos[j][2];

			for(int nx=-ncut;nx<=ncut;nx++){
				for(int ny=-ncut;ny<=ncut;ny++){
					for(int nz=-ncut;nz<=ncut;nz++){
						nq=nx*nx+ny*ny+nz*nz;
						if(nq < ncut2){  /* Spherical Sum */
							if((nq!=0) || (nq==0 &&i!=j)){
								dxn=dx+nx*lbox;
								dyn=dy+ny*lbox;
								dzn=dz+nz*lbox;
								r2=dxn*dxn+dyn*dyn+dzn*dzn;
								r=sqrt(r2);
								r3=r2*r;
								r4=r3*r;
								Co=dzn/r;
								Co2=Co*Co;
								energia+=(1.-3.*Co2)/(r2*r);

								angPart1= (5.0*Co2-1.0)/r4;
								angPart2= (5.0*Co2-3.0)/r4;
								fx=angPart1*dxn/r;
								fy=angPart1*dyn/r;
								fz=angPart2*dzn/r;

								force[i][0]+=fx;
								force[i][1]+=fy;
								force[i][2]+=fz;

								//force[j][0]-=fx;
								//force[j][1]-=fy;
								//force[j][2]-=fz;
							}
						}
					}
				}
			}
		}
	}

	energia*= 0.5;
	return energia;
}

double CoulombEnergyForcesReplication(double **pos, double *Q, double lbox,  int npart)
{
	double distancia(double **pos, double lbox, int i, int j, double *dx, double *dy, double *dz);
	double r4,r3,r2, r,dx,dy,dz, Co, Co2;
	double dxn,dyn,dzn;
	double angPart1,angPart2;
	double qij;
	double energia=0.0;
	int ncut=50;
	int ncut2=ncut*ncut;
	int nq;

	for(int i=0; i<npart; i++){
		for(int j=0; j<npart; j++){

			dx=pos[i][0]-pos[j][0];
			dy=pos[i][1]-pos[j][1];
			dz=pos[i][2]-pos[j][2];
			qij=Q[i]*Q[j];
			for(int nx=-ncut;nx<=ncut;nx++){
				for(int ny=-ncut;ny<=ncut;ny++){
					for(int nz=-ncut;nz<=ncut;nz++){
						nq=nx*nx+ny*ny+nz*nz;
						if(nq < ncut2){  /* Spherical Sum */
							if((nq!=0) || (nq==0 &&i!=j)){
								dxn=dx+nx*lbox;
								dyn=dy+ny*lbox;
								dzn=dz+nz*lbox;
								r2=dxn*dxn+dyn*dyn+dzn*dzn;
								r=sqrt(r2);
								energia+=qij/(r);
							}
						}
					}
				}
			}
		}
	}

	energia*= 0.5;
	return energia;
}
